Relativistic Structure, Stability and Gravitational Collapse of Charged Neutron 

Stars 



in 
o 
o 

(N 

» 

o 
O 
in 

(N 



o 



vn 
o 

u 



X 



Cristian R. Ghezzi 
Department of Applied Mathematics, 

Instituto de Matemdtica, 
Estatistica e Computagdo Cientifica, 
Universidade Estadual de Campinas, 
Campinas, Sao Paulo, Brazil 
(Dated: February 7, 2008) 

Charged stars have the potential of becoming charged black holes or even naked singularities. It is 
presented a set of numerical solutions of the Tolman-Oppenheimer-Volkov equations that represents 
spherical charged compact stars in hydrostatic equilibrium. The stellar models obtained are evolved 
forward in time integrating the Einstein-Maxwell field equations. It is assumed an equation of state 
of a neutron gas at zero temperature. The charge distribution is taken as been proportional to the 
rest mass density distribution. The set of solutions present an unstable branch, even with charge 
to mass ratios arbitrarily close to the extremum case. It is performed a direct check of the stability 
of the solutions under strong perturbations, and for different values of the charge to mass ratio. 
The stars that are in the stable branch oscillates and do not collapse, while models in the unstable 
branch collapse directly to form black holes. Stars with a charge greater or equal than the extreme 
value explode. When a charged star is suddenly discharged, it don't necessarily collapse to form a 
black hole. A non-linear effect that gives rise to the formation of an external shell of matter (see 
Ghezzi and Letelier 2005), is negligible in the present simulations. The results are in agreement 
with the third law of black hole thermodynamics and with the cosmic censorship conjecture. 

PACS numbers: 04.25.Dm;04.40.Nr;04.40.Dg;04.70.Bw;95.30.Sf;97.60.Bw;97.60.Lf 



I. INTRODUCTION 

The study of charged relativistic fluid balls at- 
tract the interest of researchers of different areas of 
physics and astrophysics. There exists a general con- 
sensus that astrophysical objects with large amounts 
of charge can not exist in nature m, [3. This point 
of view had been challenged by several researchers ^ , 
[27], 2^], [23. It cannot be discarded the possibil- 
ity that during the gravitational collapse or during an 
accretion process onto a compact object the matter 
acquires large amounts of electric charge. This has 
been considered in |^ and f37'|. 

In this paper we study the stellar structure and tem- 
poral evolution of compact charged fluid spheres, in- 
dependently on the mechanism by which the matter 
acquires an electric charge. 

From a pure theoretical point of view, the collapse 
of charged fluid balls and shells are connected with 
the laws of black hole thermodynamics. The third law 
of black hole thermodynamics states that no process 
can reduce the surface gravity of a black hole to zero 
in a finite advanced time 0, t33], 26]. The surface 
gravity, k, plays the role of a temperature while the 
area of the event horizon is equivalent to the entropy. 
Israel |0 gave a formulation and proof of the third 
law . It was demonstrated that the laws of black hole 



thermodynamics are analogous to the common laws 
of thermodynamics '2^. Translated to the physics of 
fluid collapse, the third law implies the impossibility of 
forming an extremal black hole, for which k = 0. An 
extremal black hole has a total charge Q : Q = VGM, 
where G is the gravitational constant, and M is the 
total mass of the black hole. The extremal black hole 
has not an horizon and constitutes a naked singularity, 
contradicting the cosmic censorship hypothesis. So, 
in this paper we explore the possibility of forming an 
extremal black hole from the collapse of a (charged) 
compact object. We will show that the answer is no. 

The relativistic equations for the collapse of a 
charged fluid ball were obtained by Bekenstein |^. 

So far in the literature, as long as the author knows, 
it was studied the dynamics of charged shells or scalar 
fields falling onto an already formed charged black 
hole (see U 0, 112, and references therein ), a nd so- 
lutions of the Tolman-Oppenheimer-Volkoff J3(l| equa- 
tions for charged stars (see 0, [s^l, ^J, etc). 
However, in the stellar collapse, the collapsing matter 
determines the background geometry and the metric 
of space-time is changing with time as the star col- 
lapses. Moreover, in general relativity, the pressure 
of the fluid contributes to the gravitational field and 
strengthen it. From the present study it is possible 
to understand if the Coulomb repulsion will prevent 
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or not the total collapse of a charged fluid ball, and 
which are the stability limits for a charged relativis- 
tic star. We emphasize that we will not be concerned 
here with the mechanism of electric charge generation, 
nor with effects that could neutralize or discharge the 
star. 

In the present paper, the equations for the evolution 
of charged fluid spheres in Schwarzchild-likc coordi- 
nates are presented in Sec. ^1 In the Sec. IIIII we re- 
derived the equations of hydrostatic equilibrium and 
the relativistic equations for the temporal evolution of 
charged spheres in a form closer to that obtained by 
May and White [l^ (for the case of zero charge). 
The equations obtained are well suited for performing 
a finite difference scheme, and allows a direct imple- 
mentation of very well known numerical techniques 
[l3], 22], 23]. The formalis m is also compatible with 
the Bckenstein equations Q , and with the Misner and 
Sharp equations 1241] for the case of zero charge. 

The calculations performed in this study are rather 
lengthy in general and we tried to be the more self- 
consistent as possible. However, we apologize that in 
most of the paper we only give hints for the compu- 
tation of the intermediate algebraic steps in order to 
keep the paper at a reasonable length. 

The matching conditions between the exterior and 
interior solutions of the Einstein-Maxwell equations 
are given in the Sec. IIVI 

For completeness, in Sec. 0it is considered an equa- 
tion of state of a zero temperature neutron gas. 

Two codes were built and are used in the present 
study: one to obtain neutron stars models in hydro- 
static equilibrium, presented in the Sec. IVII and the 
other code is to evolve the stars forward in time, in- 
tegrating the Einstein-Maxwell equations. This code 
is introduced in Sec. IVIII In all the simulations, the 
charge distribution is taken proportional to the rest 
mass distribution. 

We found that although the total binding energy 
grows with the total amount of charge, the binding 
energy per nucleon tends to zero as the charge tends 
to the extremal value. This is related to the impos- 
sibility of obtaining bounded solutions of extremely 
charged fluid configurations. In addition, it results 
that the upper limit for the total amount of charge 
could be slightly lower than the extremal value. This 
is discussed in the Sec. IVIIII 



In the Sec. IVIIII is discussed the stability and evo- 
lution of the stellar models. We study, as well, the 
evolution of stars that were suddenly discharged to 
take into account the case of charged spheres being a 
metastable state of more complex scenarios. 

In a recent work Ghezzi and Letelier ^3] presented 
a numerical study of the collapse of charged fluid 
spheres with a polytropic equation of state, and with 



an initial uniform distribution of energy density and 
charge. They found that during the evolution a shell 
of higher density is formed near the surface of the im- 
ploding star. From that work wc can not say if a shell 
will always form when changing the initial conditions 
and the scenario for the collapse. In the present study, 
we checked that the effect is negligible during the col- 



lapse of a charged neutron star (see Sec. IVIII|) . 
In the Sec. IIXI we end with some final remarks. 



II. RELATIVISTIC EQUATIONS 

Considering a spherical symmetric fluid ball, the 
line element in Schwarzchild-like, or standard form, is 
(see m, 0): 



""dt^ -f e^dr^ + r'^de'^ + r^sin^ddcf^ (1) 



An observer moving radially have a 4-velocity 
w'' = (it°, M^, 0, 0), and the electric 4-current is 
— (j'^j 0, 0). The energy-momentum tensor of 
the charged fiuid is given by: 

T"" = [5 + P)u^'u'' + P g^''' + 

^[F'''K-\9^''F'^^F^p\ (2) 

where 5 is the density of mass-energy, P is the scalar 
pressure, and F'^" is the electromagnetic tensor. The 
electromagnetic field satisfies the Maxwell equations: 



and 



(3) 



(4) 



Only the radial component is non-zero, and the 
last equation is satisfied if — — The covariant 
derivative in Eq. ^ can be written as 



(-5) 



with (-5)1/2 ^ p^{\+u)/2^ Defining 

a = i(A-|-z/) , 
the Eq. © gives 3 



dr 



and 



dt ■' 



(5) 



(6) 



(7) 



(8) 
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Integration of Eq. (TJ gives Q 

f^oi =e-"Q(t,r)/r2, 

where 

Q{t,r) = [ Airr^fe^dr, 
Jo 



and from Eq. ^ 



dQ A 2-1 a 

—— — —Anr 7 e , 
dt 



(9) 



(10) 



(11) 



From this equation we see that the charge is a constant 
e outside the fluid ball, where =0. Far away from 
the sphere, where e" — > 1 (see below), the Eq. (O 
gives the electric field of a classical charged particle. 
The Einstein's equations to be solved are: 



2^ c4 



(12) 



here (and only here), R denotes the scalar curvature, 
and R^'^ is the Ricci tensor. 

The components of the Einstein equations are (see 
i, [23): 



TO : ^ - Stt [{5 + P) u^'uo + P] 



^0 • ^4 



A' 1 



(13) 



Tl:-^- Stt [{5 + P) u^m + P] 



(14) 



: 87r((5 + P)u^uo = e-^A/r. (15) 

Multiplying the right hand side of Eq. I|13|) by 
and rearranging, the equation can be cast in the form 



-'^^^^^ + l = %-%Tr [[5 + P)u%o + ■ (16) 
dr 

Setting e^^ = 1 + f + g, with / = — 2m/r and g — 



(17) 



the Eq. becomes 



(18) 



Integrating this equation, we get the equation for the 
mass 0: 



47r 



C ./O 



[((5 + P) u^ua + P] dr 

1 1 



2c2 r 2c2 



dr + ma, (19) 



where mo is an integration constant, and can be taken 
as zero for the purpose of this work. 

Outside the fluid ball, the mass do not depend on 
r and then: m : m{t,rs) for r > r^, where is the 
coordinate r of the surface of the sphere. Actually the 
mass do not depend on t for r > r,; from Eqs. (|15|l 
and (fT7|l : 

^^9.§+,,(^S + P)r^uou\ (20) 
dt r dt 

as Q is constant and 5 = _P = outside the fluid ball, 
TO is independent of t. 

Subtracting Eq. (fT!^ from Eq. (fTHl it is 

— ^^^-^=87r((5 + P)(7.iui-u%o). (21) 
r ctr 

So, A + 1/ is independent of r. To get an asymptotic 
flat solution it must beA + i^^Oasr^oo. Thus a 
solution of Eq. if^ is 3] 



A = —V for r > 



(22) 



Using the Eq. ifTzll . if^. and substituting into Eq. 
(Q, the line element outside the ball is 



ds^ = -{l-^+QL-)dt^ + {l-^-^ + %)-^ dr^ 
r"sm-fd0', (23) 



which is the Reissner-Nordstroni space-time. The 
gravitational mass is M = m{rs) = constant. We 
see that although the (interior) mass distribution and 
the metric depends on time, the external space-time 
is static. This is the Birkhoff theorem, and implies 
that a spherical distribution of mass and charge can- 
not emit gravitational waves. 



III. EQUATIONS IN CO-MOVING 
COORDINATES 



In this section we will specialize the equations for 
the case of a coordinate system co-moving with the 
fluid. We will use a notation closer to that used in the 
May & White papers p^ . 

The 4-velocity for an observer co-moving with the 
fluid is = (a~^, 0, 0, 0), and the measured 4-current 
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is = (j'^, 0,0,0). The 4-velocity satisfies: u'^u^ = 

The line element for co-moving coordinates in 
"standard" form (see Pag- 336) is: 

ds'^ = a{t, ^ifc^df - 6(t, ^xfd^'^ 

~R{t,^if{de + iiin^ed(j)'^) , (24) 

note that now a, 6, and R are functions of /i and t. 
We will see in the next subsection that the coordinate 
fi can be chosen to be the total rest mass inside a 
sphere, such that each spherical layer of matter can 
be labeled by the rest mass it contains. 

The stress tensor in co-moving coordinates is ob- 
tained using the solution of the Maxwell equations 
®: 





= 5^ + 


87ri?4 


(25) 


n 




87ri?4' 


(26) 


T| 


= r| = 


87ri?4 


(27) 




- - 


0. 


(28) 



From now on we will restore in the equations the grav- 
itational constant G and the speed of light c. 

Similarly as was defined above, 5 (? is the density 
of energy in [dyn/cm^] and P is the scalar pressure in 
the same units. It will be useful to perform the split 

m 

5 = p(l + e/c2), (29) 

where p is the rest mass density and p e is the internal 
energy of the gas. 

A. Equation of particle conservation 

We want to set the radial coordinate equal to the 
rest mass enclosed by co-moving spherical layers. We 
performed a calculation analogous to that of May and 
White 1 and used the same notation for clarity and 
comparison. 

In this work is assumed that there is only one 
species of particles, and there are no particles created 
nor destroyed. The number density of particles will 
be denoted by n and the rest mass density by p. If 
the rest mass of the particles is m„ then: 

p — n rrin . (30) 

The conservation of the baryon number current is ex- 
pressed : 

nu'f^ = (31) 



or equivalently: 

pu';, = 0. (32) 

The line element in spherical coordinates written in 
"standard" form 43] is : 

ds^ = a{t, Rf dt^ - b{t, Rf dR^ 

-R^{de + sm^ed(l)'^), (33) 

observe that this line element is not identical with that 
given in Eq. (|^ . 

The conservation of baryons in this frame is reduced 
to (see Eq. E21): 

(pi?2&)^, = 0, (34) 
thus the quantity {p R^ b) is a function of R only 

{pR'b)^f{R). (35) 
Choosing f{R) = l/4 7r: 

b=l/4TTpR^. (36) 
On the other hand, the proper mass is 

p = [ p^gd^x, (37) 

JV 

where \/ ^^^gd^x — b R^ sin 9 dR d9 d(f> is the volume 
element of the space section , and V is the volume 
of integration. Then 

H= AnpR^bdR, (38) 
Jo 

is the proper mass enclosed by a sphere of circum- 
ference 2tt Rs, and radial coordinate Rg ■ With this 
election of coordinates: 

dn = dR. (39) 

So we can perform the coordinate transformation 

(t, i?, 0, 0) (i, /i, 0, 0) , 

and in this case the metric functions change to 

a{t,R) ^ a'{t,p) , (40a) 
b{t,R) ^b'{t,fj.), (40b) 
R^R'{t,p). (40c) 

Replacing these functions on the line element H33(l . 
dropping the primes, and using the Eq. 139|) we obtain 
the line element given in the Eq. (|24|l . So, the radial 
Lagrangian coordinate was gauged to be the rest mass 
"/i" of each spherical layer of matter. 
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B. Charge conservation 



The proper time derivative of the charge is (see Eq. 



where m is given by the Eq. H45|) . Combining Eqs 
(|36|) and (|42|) we get the mass conservation equation 



(49) 



= 0. 



(41) 



because it is possible to write the 4-current as a prod- 
uct of a scalar charge density times the 4- velocity, i.e.: 
j'^ = pch (m°, 0,0,0). Thus the electric charge is con- 
served in spherical layers co-moving with the fluid. 



C. Einstein-Maxwell equations in the co-moving 
frame 



The Einstein equations (see the Appendix) for the 



It can be seen that the quantity u (see Eq. I47II is 
the radial component of the 4-velocity of the fluid in 
Schwarzschild coordinates. Using the transformation 
rule for the contra- variant 4- vector u"^ from co-moving 
to Schwarzschild coordinates: u'' 



{dx"'/dx'^)u" 



yields 123: 



cTt/a 



Rt/a. 



components G* = are 



—R.t + -T'R p- ~ ^.A't = ^ ' 
a 



(42) 



where we use the notation R t ~ dR/dt, R^^t = 
d^R/dtdp, etc. 

The equation for the component G* is (4^: 



D. Conservation of the energy and momentum 

From the Einstein equations and the contracted 
Bianchi identities follows that: 



AnGSR^Rf, = — 
2 

GQQ,^ 



R 



RRt ■ 



RR,n ^ 
62 



GQ2 



c'^R 



(43) 



In co-moving coordinates the components are: 



^ -.11. 



b 



R 



P. 



and the equation for the component G^^ is|46 



((5c2- 



P) 



47rG 



-PR^Rt = - — 



R 



RRl RR\ GQ^ 



62 



c^R 



Am ((5c2 -h P) i?4 



= 0. 



(50) 



0,(51) 



(52) 



= -Gmt . 



(44) 



Expanding the Eq. and using the Eq. l|5T|) we 



get: 



We introduced into Eqs. (|^ and (|^ . the defini- 
tion of the total mass: 



(53) 



1 OO 

m{p,t)^ATi I SR'^R,f,dp + — -^^d^. 
^0 ' c Jq R 



(45) 



We give further definitions that simplify the aspect 
of the equations: 



r = ^ = AnpR'R^,, 
R.t 



(46) 
(47) 



which is identical to the non-relativistic adiabatic en- 
ergy conservation equation, and express the first law 
of thermodynamics. It could be surprising for the 
reader that there is not an electromagnetic term on 
this equation: this is due to the spherical symmetry 
of the problem. 

The Eq. can be written: 



With this definitions an integral of the Eq. H43|l and 
(gH) is: 



[aiS + P/c')/p]^^ 
aiS + P/c^)/p 



1 



[(5 + P/c2)/p]c2 



P 



AttR^P 



(54) 



= 1 + ^ - 



2mG Gg2 



i?c2 c4i?2 ' 



(48) 

or using the Eq. and the definition of the rel- 
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ativistic specific enthalpy (see [l^, |23|) w — {S + 

P)/p^: 



e P 
w = 1 + — ^ 



(55) 



we can write the Eq. (|54|l in a form closer to that used 
by May & White ;22j: 



1 



aw wc-^ 



(56) 



E. Equation of motion 

The equation of motion can be obtained expanding 
the Eq. ll^ . using the Eq. (g^l to ehminate the 
factor R.^ti the Eq. H52|l to eliminate a^/a and the 
definitions given in the Eqs. 
After lengthy calculations we obtain: 



u.t 



-a AnR^ 



r 



p, 



C2 i?3 



Gm 



(57) 



This equation reduce to the equation of motion ob- 
tained by Misner and Sharp and May and White 
[2^ for the case Q — and is compatible with the 
equation derived by Bekenstein JJ. The Eq. H57() re- 
duce to the Newtonian equation of motion of a charged 
fluid sphere letting c — > cxd, and consequently: F — > 1 
(see Eq. 021; w ^ 1 (see Eq. EHJ; and a — > 1 (see Eq. 



F. Equation of Hydrostatic Equilibrium 

The equation for charged fluid spheres in hy- 
drostatic equilibrium is a generalization of the 
Tolman-Oppenheimer-Volkov (3Qj equation obtained 
by Bekenstein Q . We can obtain it from the equation 
of motion (|57|) . taking m = 0, and u^t = 0. In hydro- 
static equilibrium, the factor F (see Eq. 1^ is 



F^ = 1 - 



2mG GQ^ 



Rc^ d^R^ 



(58) 



Rearranging 47J the Eq. (|57|l and using the 
Eqs. H55|l and (|58|) . the Tolman-Oppenheimer-Volkov 
(TOV) equation for charged fluid spheres is: 



dP 2 
dR ^ 

Q dQ 



'niG 



4ttG 



P) 



PR^ 



RlR- 



2mG 



47ri?4 dR ' 



(59) 



This equation gives the well known Tolman- 
Oppenheimer-Volkov "scJl equation when Q = 0. 



IV. THE MATCHING OF THE INTERIOR 
WITH THE EXTERIOR SOLUTION 

In the Sectional we find an exterior solution to the 
problem which agrees with the Reissner-Nordstrom 
solution. In the preceding section we derived an in- 
terior solution that we want to integrate with a com- 
puter. However, we must show that the interior and 
the exterior solutions match smoothly. The matching 
conditions between the interior and the exterior solu- 
tions are obtained by the continuity of the metric and 
its derivatives along an arbitrary hypersurface. First 
of all we must establish the equality between the exte- 
rior and the interior metric at the surface of the star. 
This can be done transforming the interior solution 
given in co-moving coordinates into the Schwarzschild 
frame. The metric tensors in the two frames are re- 
lated by the tensor transformation law: 



(60) 



The line element in the Schwarzschild frame is given 

by 

ds^ = dT^ - dR^ - i?2 dVi^ , (61) 

with d^ = dO + 8\ii^9d4>. Using the Eq. we ob- 
tain the non-trivial relations between the metric coef- 
ficients: 



1 



1 

~B2 



T2 ( 



1 







R.t T.t 



R 



- R.u T 



(62a) 
(62b) 
. (62c) 



From Eq. I|62b|l and using Eqs. (gSl and (|T7|l we get: 



or equivalently. 



B^ 1 



2mG GQ^ 



2\ -1/2 



i?c2 C^R^ 
From Eqs. and (|^ : 

1 fTiVfT^-u" 



(63) 



(64) 



(65) 
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With a little amount of algebra we can transform this 
equation into: T"]. = B'^V'^a?, and using the Eq. 



(66) 



This last equation can be obtained independently by 
comparing the induced metric, from interior and exte- 
rior solutions, on the hypersurface S: t = t, fi = — 
constant, 9 = 6, 4> = (f). The line element correspond- 
ing to the exterior solution is given by: 



dsl = dT^ - B^ dR^ - R\ dn\ , 



(67) 



from now on a plus (minus) subscript or superscript 
denotes exterior (interior) solutions. The line element 
for the interior solution is (see Eq. I24|l : 



Due to the spherical symmetry: dO- 



ds^_ = dt^ - dn"^ - Rl dnl . (68) 

dn+. The 

line element compatible with the induced metric on 
the hypersurface, from the exterior solution isfisf: 

ds^lj. = {A^ rj - B^ i?2 ) dt^ - R\ dnl , (69) 

and the metric induced on E by the interior solution 
gives: 

dst^^ ^ dt^ - R^_ dni . (70) 

On S we have rfs^|j, = ds^|j^, so we obtain: 

R+ = R- (71a) 
A^ (? - i?2 ^a^(? . (71b) 

and we recovered the Eq. (jSHl)- So Eq. (|71b(l (or 
equivalently Eq. I66|l are conditions for the continu- 
ity of the metric along E. From now on, we will use 
indistinctly R — R^ = R~ . The relation between 
the exterior and interior metric coefficients constitutes 
three equations with the four unknowns (see Eq. I62|) : 
R,t; R,fi; T^t', T^^. Observe that R t and i?.^ are ob- 
tained as part of the interior solution. In addition, 
the solution exterior to the sphere of matter must sat- 
isfy the field equations: R^u = 0, here is the Ricci 
tensor. So we obtain (see for example ,3^ Pag. 180 
and Pag. 337): 



A{R) c - 



B{R) 



(72) 



The continuity of the metric is almost already estab- 
lished. Only for completeness, and after a little alge- 
bra, we can find an equation for T^: — h"^ B^ 
(although we will not use it). 



The equation of the hypersurface containing the 
surface of the star is 



: t ^ t, ^ = fis 



</), (73) 



where /is is the Lagrangian mass at the surface of the 
star, so the metric coefficient B is: 



B= 1 



2MG GQ"' 



2\ -1/2 



(74) 



with M = m{t,fis), R = R{t,iis), and Q — Q(t,^s)- 
The coefficient A^c^ = gu is A'^c^ — 1/B^, and we 
recover the Reissner-Nordstrom solution: 



1 

dR^ - R^d0^ 



1 - 



2MG GQ^^'^ 



Rc^ c'^R^ 

2 



sin^ 



It is important to remark that we obtained the 
Reissner-Nordstrom solution by performing a tensor 
transformation of the interior solution. The exterior 
solution must satisfy the Einstein field equations, and 
so we used them to find the Eq. (|72|l . Thus, the 
exterior and interior metric are consistently matched. 

It remains to prove the continuity of the derivatives 
of the metric along Ss . This is a more lengthy calcula- 
tion and we will give only the the main algebraic steps 
and the final results. The continuity of the derivatives 
of the metric is established from the continuity of the 
second fundamental form, or extrinsic curvature ten- 
sor Kab of the hypersurface By definition (see 
Pag. 59): 



a P 
T'a;[3 ^a^b ' 



(76) 



here n" is a unit vector normal to S^; e" = dx"" /dy°', 
are basis vectors on Sg; the coordinates on are 
denoted with while x°' are the space-time coordi- 
nates. Equivalently : Kab — \ {^ngap)e.ae^, here £„ 
is a Lie derivative. 

We will choose = {t,0,(j)) as coordinates on E^. 
The equation of the hypersurface approaching from 
the exterior is: 

S+ : T = Tit, Ais), i? - R{t, ^s), ^ - ^, <^ = <^ , (77) 

The basis vectors on are: 



ef,) =a-i(Tt,i?,t,0,0) 
e?g.^{0,0,R'\0) 



= (0,0,0,i?"isin~^6') , 



(78a) 
(78b) 
(78c) 
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in the coordinate basis B = {dT,dii,d0,d^}, and the 
normal co- vector is: 



{-Rt,Tt,Q,0). 



(79) 



The equation for Es as seen from the interior is (see 
Eq. CSll: 



: t = t, fi = ^s, 9 = 9, (j) = 
the basis vectors on are: 

„a _ /„-l 



e^,) = (a-\ 0,0,0) 
e^,) = (0,0,i?-\0) 



e^^) = (0,0,0,i?~W 



(80) 

(81a) 
(81b) 
(81c) 



in the coordinate basis B' = {dt,df_i,de,dcj,}, and the 
normal co- vector is: 



n- = (0,6,0,0) 



(82) 



If there are no surface distributions of energy-matter 
on Es, it must be 49]: 



(83) 



The angular components of this equation are Kgg = 
Kgg, and = K^j;^. For K+ we have: 

Ke - ^eenB^efg^elg) (84) 

here F^^ is the connection. From now on, in the 
present section, we are not using the sum rule over 
repeated indexes. In the last two steps we used the 
Eqs. f7TH|) and (|72J). 

For Kgg, wc find Kgg = R R^^ /b, and using the 
Eqs. (gg and IHI): Kgg = F", so Eq. ^ gives: 



F+(i,/i,,) = F-(i,M.): 



(85) 



or equivalently: 



1 + ^- 



2MG ^GQ"" 
~R^ ^ 

2ms G 



i?2 



-1/2 



Rc? 



C4 i?2 



2\ -1/2 



(86) 



Here M = m(i, /ig -f e) and Q = (5(t, + e), with e a 
positive arbitrary small number or zero. The defini- 
tion of u assures its continuity across E^, so u"*" = u~ . 
Moreover, we already found that M and Q are con- 
stants independent of t and /i, so the Eq. 186(1 implies 



that rhs = 0. Then, using the Eq. ((20|l in co-movcl co- 
ordinates (Q = in both coordinate systems, see Eq. 
inl and Eg. I41II . we obtain the boundary condition at 
the surface [sol: 



Ps = Pit, ^is 



0. 



The equation [if^^] — gives an analogous result. 
The component K^^ is: 



tt -2 

a 



a b 



(87) 



The external component is: 

J- - I t t t T 

J^tt ~ "-t:t '^(t) "^(t) + "-t;!- 6(4) e(f) + 

rir-t e\t) ejt) + n^-r 6(4) e^j^ , 

the sum rule is not used here. For clarity, we will 
define <? = = go Eq. I|71b|) becomes 



(89) 



and an over-dot means dt- Applying the over-dot op- 
erator to this equation we obtain an equation that will 
be used later: 



So, 



2f-'^RR-2f^f + 2aa + a^f-'^f 
2t7 



(90) 



K+ - {- 
{- 
{ + 
{+ 



d{a-^R) I ,df . _i 



tt L QJ. 

djq-^R) 
OR 

dT 

dR 



lf-'^ia-'R)}a-'fR + 
\r'^{a-'R)}a-'fR + 



2^ dR 



{a-^f)}a-^R^ 



(91) 



Expanding this equation, and using the Eq. I|90() to 
simplify the result we obtain: 



K^,=a 



-3 



2aaR^ - 2a^ RR- a'^ f 
2tfR 

further simplifications can be made using: 



• 2mG • 2Q2G 2mG 

J = „o o „^ . R- 



i?2 c2 i?3 c4 

ii Rd 
/T = a(i?2/a2 + /)i/2 = ar. 



U = 



(92) 



(93) 

(94) 
(95) 
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Thus we obtain: 



mG GQ^ 



Gm 



R^c^ c^R^ RRc^ 



(96) 



Comparing this equation with the Eq. (|87|l we obtain 
an expresion for [Ku] ~ 0: 



u 
a 



iG GQ^ 



Gi 



i?2 



c^R^ RR 



a b 



(97) 



Now rearranging and remembering that 5 = 1/4 tt p R^ 
and Gm/R = -Air G P R"^ /c^ , using the Eq. ^ to 
replace a^c^/a we get: 



Ut 

AttG 



i?2 



PR 



GQ^ 

C2 i?3 



(98) 



Then, we see that the equation [Ku] = is satisfied 
by virtue of the equation of motion (see Eq. I57|l . and 
thus the exterior solution matches with the interior 
solution. 



V. EQUATION OF STATE 

The equations obtained above must be supple- 
mented with an equation of state for the gas (EOS). In 
particular we choose to represent the gas of neutrons 
as obeying the quantum statistic of a Fermion gas at 
zero temperature. This let us use the Oppenheimer- 
Volkov numerical results as a test-bed for the code. 
However, any other EOS can be equally implemented 
in the two codes presented in this paper. 

The energy per unit mass is given by 



1 dx . 



(99) 



where x = p/m„c is the relativity parameter, and p is 
the momentum of the particles. 

The Fermi parameter is: xp = Pf/"t.„c, with pp = 
^^n, where n is the number density of neutrons. It 
can be written (we suppress the sub-index F from now 
on): 



1/3 



(100) 



where 



Po 



37r2ft3 



= 6.10656 X 10^^ [gcm-3]. (loi) 



The Eq. H99|l can be integrated to give: 
m^c^ 1 



5c' ^ 



87r2 



■ \/TTx2 (l + 2^2) 



log (a; -I- i/l-f a;2) 



(102) 



the quantity S is measured in [ergcm"^]^ or equiv- 
alently in [dyncm~2j. This is an expression for the 
total energy of the gas, so it includes the rest mass 
energy density. Although this is obvious by definition 
(see Eq. I99|) . it can be useful to check this expanding 
Eq. (|102|l for a; ^ (or x ^ oo). In this two cases an 
split similar to Eq. (|29|l is obtained. 
The pressure is given by ^36; ]: 



P 



3^2 ft3 

and performing the integral 

4 5 1 

r7i„ c 1 



dx , 



(103) 



P = 



ft3 87r2 



Vl + a;2 (2xV3- 1) + 



log ((a; + vT+a?2) 



(104) 



measured in [dyncm~^]. 

For charged matter, the equation of state must de- 
scribe several species of particles, i.e., neutrons, pro- 
tons and electrons. However, the amount of charge 
in all the models studied here is very small to change 
the EOS. For an extremal charged fluid ball there are 
only one unpaired (not screened) proton over 10^^ par- 
ticles. In this case, the change on the chemical poten- 
tial of the neutrons is negligible and charge neutr ality 
can be assumed from a microscopic point of view |5ll | . 

The simulations were extended beyond the densities 
at which the EOS is valid. However, we want to take 
into account possible charge and pressure regeneration 
effects, that occurs at high densities. 



VI. NUMERIC INTEGRATION OF THE TOV 
EQUATIONS FOR CHARGED STARS. 

The equations that must be integrated to obtain 
the stars in hydrostatic equilibrium are conveniently 
written in dimensionless form. From Eqs. |^ and 
(|16|) it is possible to obtain an equation for the metric 



component — b~ 



dX 
dR 



GQ^ 

c'^R^ 



SttG 



R5- 



R 



(105) 
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In addition we use the chain rule on the TOV equa- 
tion (Eq. EHJ: 



since po is the factor in the expressions for the pres- 
sure (Eq. I103|) and energy (Eq. I99|) . Hence: 



dp 
dR 



(niG I 47rG 



-((5c2+P)- 
Q dQ 



PP3 „ G<3 



1 



1 



47ri?4 rfi? {dP/dp) 



(106) 



The dimensionless equations (|106|) and IjlOSI) can be 
obtained with the replacements: 

m = morn, P = PqP, Q = QoQ , P ^ P-o P' 

R = RqR, 5 = SqS , Pch — Pcha Pch ■ 

Where the bar denotes dimensionless variables, and 
the subscript "0" indicates the dimensional constants. 

The full set of dimensionless equations to be inte- 
grated are: 



dX_Q^ 
dR~R^ 
dQ 



e^R5 



1 



R 



dR 



p. 



•CHC'/'R' 



^^pe'/'R' 
dR ' 



(107a) 
(107b) 
(107c) 



dp 
dR 



---i'S + P) 

Q dQ 1 
R^ dRdP/dp 

^^SR'dR+S-^. 
dR i?2 ^ji 



(m + _ ^) ^ 
R{R ~2m+^) dP/dp 



(107d) 
(107e) 



The dimensional constants are given by the set of 
equations: 



Afo =47rpo R'o 

Po = Mq 

Qo = VGMo 
„ MoG 

-Ko — 5- 

Po = Po 

Pcha Po 
So = Po- 



V4 7rpo G3 



(108a) 

(108b) 
(108c) 

(108d) 

(108e) 

(108f) 
(108g) 



Mo = 1.03706 X 10^3 g 
Ro = 7.69944 x 10^ cm 
p^^^ = 1.80808 X 10" gcm"^' 
poc^ = 1.62502 X 10^^ dyncm"^ 
Qo = 2.67888 x 10^^ StatCoulombs . (110c) 



(110a) 
(110b) 
(110c) 
(llOd) 



Moreover the set of equations (|107a|) - (|107e|) are 
invariant under the transformations: Rq — > Rool, 
Mq — > Mga, and po Po/a^j etc (the other con- 
stants can be obtained from these), with a an arbi- 
trary number different from zero. So, it is possible 
to choose any other set of constants consistent with 
these transformations. In particular, the dimensional 
constants obtained by Oppenheimer and Volkov |30| . 
are recovered by setting: a = v327r2. jj^ ^j^jg case: 
Ro = 1.36831 X 10^ cm; Mq = 1.84302 x lO^"* g, 
e tc, equivalent to the constants obtained in that paper 



In the Eq. (|107d|l the term dP/dR, is obtained 
deriving the Eq. H103|l : 



dP 
dR 



37r2x/r 



(111) 



It is used a 4*'' order Runge-Kutta scheme to in- 
tegrate simultaneously the set of equations H107a|l - 

The charge distribution is chosen proportional to 
the rest mass distribution: pch = o:{p,t) p. For sim- 
plicity, in this paper we will concentrate on the case 
a = constant. 

The code implemented to integrate the equations 
above is called HE05vl. The IIE05vl is used to build 
neutron star models in hydrostatic equilibrium. This 
models constitutes an initial data set for another code 
presented below. 



VII. NUMERICAL INTEGRATION OF THE 
EINTEIN-MAXWELL EQUATIONS 



It is natural to choose 

4 3 



X 10" gcm-3 , 



(109) 



In this section we collect the equations that must be 
integrated numerically in order to simulate the tem- 
poral evolution of the stellar models constructed with 
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the HE05vl code. The set of equations is: 



du = 



4 7ri?^ 



P., 



QQ., 



AttG 



PR^ 



w \ ■ inR-^ 



c2 i?3 



dt 



w = 1 + ^ 



2 2toG GQ2 



aw = aoioo exp 



de + Pd 



dr7i = 47rp(l + -j)i?2 R^^ dp 



1 

c2 i? 



pi? = pqR^ exp 







dg = 4 7rpch R^dR/T 
dp = A'KpR^ dR/T . 



G m 
~R2~ + 

EH (112a) 
EH (112b) 

EHl (112c) 
|ll](112d) 

f 

ISni (112e) 

m (ii2f) 

ISl(112g) 

|SSl(112h) 

(112i) 

m (ii2j) 



Where we mean oqWo = (aw)(i, 0) and PoR^ ~ 
(P i?^) (0, ^) . The references to other equations are en- 
closed in square brackets. Comparing the Eqs. (|112|l 
with the May and White equations (231, it can be 
seen that they share a very similar structure when 
Q = 0, y{p,t). We have chosen a charge distribu- 
tion proportional to the rest mass distribution so we 
obtain the Eq. ifTT^ . The Eq. ((11^ is the only 
equation that have a second time derivative and con- 
stitutes the dynamical part of the Einstein-Maxwell 
equations. The Eqs. Ijll2t|l and pi2g| ) represent the 
Hamiltonian and momentum constraint equations, re- 
spectively [S^l- The numerical code to integrate the 
equations above are called Collapsc05v3. 

The initial conditions at t = to are obtained with 
the HEOSvl for each stellar model, and consist on: 
the initial mass density distribution p{p,tQ); the ini- 
tial electric charge density distribution Pch{f-i ^o); the 
kinetic energy e(/i,to); the cell spacing dR{p,tQ) and 
the surface radius R{ps, to)- Here, ps is the mass coor- 
dinate of the surface. Then, the output of the HE05vl 
is used as an initial Cauchy data set in the space-time 
hypersurface t — t^, p = p, 9 ~ 9, (p ^ (p, for the 
Collapse05v3. 



The boundary conditions are: 

P = 0, eX p = ps, Vt 

a = 1, at p = Ps, Vt 

u = 0, aX t ~ to, y p 

p = 0, at /i = 0, Vt. 



(113a) 
(113b) 
(113c) 
(113d) 



The condition given in Eq. I)113bp . makes the time 
coordinate synchronized with an observer co-moving 
with the surface of the star. The conditions expressed 
in the Eqs. H113a|) - (|113dp . also gives: 

r = l, Q = 0, m = 0, at/i = 0, Vt 



The finite difference algorithm was implemented 
using a method similar to that developed by May 
and White [13 • The Eqs. ((Tna|) -( [TT2jl ) were inte- 
grated with a leap-frog finite difference method plus a 
predictor-corrector step. We call this numerical code 
as Collapse05v3. The integration is iterated accord- 
ing to a desired error control. The method is second 
order accurate in space and time "s^. Each experi- 
ment was repeated with different number of particles 
and with different values of the numerical viscosity 
parameter, in order to check the convergence of the 
results. In general, good results are obtained using 
~ 200 points along the coordinate p. The compati- 
bility between the two codes is an indirect check of 
the convergence properties of the Collapse05v3, since 
the HE05vl is 4t/i order accurate. For example, at 
t = to there exists differences in the integrated mass 
between the two codes of the order of ~ 0.01 M©, 
using 200 points in the Collapse05v3. This differ- 
ence can be made < 0.001 Mq by taking, instead, 800 
points, and so on. For the integration in the time 
dimension is needed to take initially small time-steps 
dt — 10^^ s, but the code is capable of self-adjust the 
time-step so as to take small changes in the physical 
variables. Therefore, with this algorithm it is indif- 
ferent to take an initial time-step of 10"'^ s or 10~^ s, 
the code always self-adjust the time-step and prove to 
converge always to the same solution. The shocks are 
treated by adding an artificial viscosity, which is non- 
zero only on discontinuities. Its effects are to smear 
out the discontinuities over several cells, and to re- 
duce the post-shock oscillations and the numerical er- 
rors related to the leap-frog method. We experienced 
with distinct functional forms for the viscosity term, 
probing the method of May and White to be the 
best in the present algorithm. However, there are not 
strong shocks (like in Ref. 14]) in the set of simu- 
lations studied in this paper and with the particular 
initial conditions chosen. 

In the Collapse05v3 the Eqs. ^(TT^i and |[TT^ are 
not integrated in time, and the rest mass and charge 
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are constrained to be constant in the layers of matter 
(as it must be). We check that this kind of algorithm 
reduce the numerical errors. 

As was showed in the Sec. the total mass- 

energy at the surface of the star is also a constant 
m{t,iJLaup) = M = constant. The numerical errors 
have an impact on the constancy of the total mass- 
energy M, and this is the physical quantity with which 
is checked the accuracy of the simulation. In all the 
runs M is very well conserved, with the exception at 
the time when a trapped surface forms. In that case, 
very large gradients forms and serious numerical in- 
stabilities appear. Although, it is possible to manage 
the code to keep the error in mass-energy conservation 
below ~ 5 % at the time the apparent horizon forms. 
Of course, once the trapped surface forms nothing can 
change the fate of the matter inside the apparent hori- 
zon (see the Sec. I VIII B l|l . The simulation must be 
discontinued soon after the matter cross the gravita- 
tional radius because the time-step attains very small 
values. In general, when dt is roughly dt < 10^^° s, it 
is needed a very large number of time-steps to obtain 
further progress and this consume too much computer 
time. The reason is that the large gradients produce 
huge changes in the physical variables and a very small 
time-step is required to keep the errors at a low level. 
Thus, it is not possible to follow the dynamical evo- 
lution until a stationary regime is reached (or the for- 
mation of an event horizon, if ever possible). When 
the outer trapped surface forms, some layers of matter 
external to it can be collapsing or expanding and is 
not possible in general to say which will be their fate. 
However, the matter that falls in the trapped surface 
can not escape from it and this is enough to say that 
the star (or some part of it) collapsed. 

There are two different sources of perturbation in- 
troduced in the models: the first source of perturba- 
tion is of numerical origin and cames from the map- 
ping of the equilibrium models in the evolution code. 
This mapping is not exact, and there is a difference 
in the integrated rest mass and total mass between 
the two codes that is taken < 0.01 M0. This dif- 
ference can be made as small as wanted by taking 
a higher number of points to perform the integration 
with the Collapse05v3; the second source of perturba- 
tion, is introduced by multiplying the kinetic energy 
by a small constant grater than one, i.e., e — > ae. In 
the present simulations we choose a to take the values 
a - 1.01 - 1.04. 

The structure of the space-time of a collapsed 
charged star is very complex, with the possible exis- 
tence of time-like singularities and tunnels connecting 
several disjoint asymptotically flat space-time regions. 
This tunnels and the true singularities, pertains to 
the analytically extended portions of the space-time. 



and the present code is not prepared to reach that re- 
gions of the total manifold. In fact, with the present 
code is not possible to reach the Cauch y ho rizon where 
important physical effects must occur , ■ How- 
ever, from an astrophysical point of view is interesting 
to know what happens outside the apparent horizon, 
where the astronomers live. 

We will not describe the numerical methods in more 
detail because is not the objective of this paper. More- 
over, the numerical methods are very well known and 
better explained in textbooks and papers (see for ex- 
ample Ref. [III). The Collapse05v3 is an extended 
version of a preliminary code developed by Ghezzi 
(2003), earlier test-beds and results obtained with this 
code were published in collaboration 0, 



VIII. RESULTS 

A. Charged neutron stars in hydrostatic 
equihbrium 

The Fig. |^) shows the mass-radius relation for 
neutron stars with zero charge. In this curve differ- 
ent points correspond to stars with different central 
densities and total number of nucleons. The results 
are in a gree ment with the results of Oppenheimer and 
Volkov I33. 

The Fig. ^) shows the mass-radius relations for 
models with different values of the charge to mass ra- 
tio Q/VCfi. 

In the Figs. the central density of the models 
increase from right to left, and following counterclock- 
wise sense inside the spiral. 

The factor Q/VCfi is a constant in each curve of 
the Fig. (Pd), but Q/VGM varies along them. This is 
because the binding energy (/i — M)c^ is not constant 
along a curve with constant Q/VCfi. In other words, 
the binding energy varies with the central density of 
the model. This can be seen in the Fig. UK); which 
in addition shows that the binding energy vary as a 
function of the charge, as well. The total binding en- 
ergy of the stars increase with its charge (see Sub-sec. 
IVIII A2|l . For example, the maximum mass model 
with Q/VGh = 0.97, has a binding energy larger than 
the maximum mass model with Q/VCfj, = 0.8. More- 
over, in each curve, the maximum is attained at lower 
densities for higher total charge (see Fig. [2t)- 

In the Fig. (|2h) can be seen that for higher values of 
the total charge (higher Qj^fGp^ the maximum mass 
and the radius of the models became larger. 

The numbered circles in the Figs. ^ and in- 
dicates some stars that have the same central density. 
This models were evolved with the Collapse05v3 for 
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the cases Q/^/G^i = 0, 0.5, 0.8 (see Sec. IVIIIB|I . 

The Fig. ||2Jl shows the mass of the models as a func- 
tion of its central density. We found that for any value 
of the charge below the extremal value {QI\fG\x < 1), 
there is still unstable and stable branches in the so- 
lutions. As it is well known from the first order 
perturbation theory the regions of the curve where 
dM/dp < represent unstable solutions, and where 
dM/dp > there are stable or marginally stable stars 
(see [3|, HH). 

There are models that have negative binding energy 
(see Fig. Et), but they are on the unstable branch of 
the solutions (see Fig. 

The maximum mass models are indicated in the 
Fig. with a vertical bar, while the tag 1st or 

2nd show the position of the first or the second mass 
maximum, respectively. The second maximum in the 
curves, are due to a pure relativistic effect. The reason 
is that "energy has weight", paraphrasing Zeldovich 
and Novikov (41J : as the number of baryons increase 
the total mass of the stars also increases attaining the 
first mass maximum. Passing the first maximum, the 
mass begins to decrease as the pressure became softer 
at relativistic energies. However, with a further in- 
crease of the density the contribution to the "weight" 
due to the kinetic energy of the particles is more im- 
portant respect to the rest mass, and the second max- 
imum appears. 

From the Fig. ijHi), it is evident that this effect also 
takes place for charged neutron stars (see also tables 

Ulnl . 

The Fig. ^ shows the coefficient grr of the metric 
as a function of the lagrangian coordinate, for model 
number 13, for three different values of the charge to 
mass ratio Q/VG^ = 0, 0.5, 0.8. 

In the Table m are shown the results of the numeric 
integration of neutron star models without charge. 
The results are in agreement with the Oppenheimer 
and Volkov calculations. The first and second mass 
maximums are indicated on the table. 

In the tablesiniand lllll are the results of the integra- 
tion of neutron star models, with QI\fG[i = 0.8 and 
Q/\fGp, = 0.97 respectively. The Tables Ull and Iml 
let us make a quantitative comparison of the charged 
neutron stars models with the properties of the neu- 
tron stars with zero charge, given in Tabled 

In particular, the maximum mass for models with 
Q/VGh = 0.97 is 8.734 M©, and the radius of this star 
is 75 km (see Table 1111}) . For comparison, this case cor- 
responds to models with a parameter / = 0.00111592, 
in the notation of Ray et al. i35]. But they used a 
charge distribution proportional to the mass-energy 
density and a different equation of state. However, 
the results presented in this section are in qualitative 
agreement with their results. 



We must observe that to study the possibility of 
making extremal Reissner-Nordstrom black holes it 
doesn't matter if the charge distribution is propor- 
tional to the rest mass density or total energy density, 
since the binding energy per nucleon tends to zero at 
the extremal case (see sub-sec. IVIII A2|l . 



1. Extremal case 

With the HE05vl it is possible to reach the sector 
QI\fG[i > 1. Black holes with this charge to mass 
ratio constitute naked singularities. 

We found that approaching the value Q/VGh = 1, 
the mass and radius of the models tends to infinity, 
within the computer capacity. Of course, it is impos- 
sible to show plots of this results, but a new technique 
was developed that will let us study this "solutions" 
and will be presented in a separated paper Ts*! . 

The Newtonian Chandrasekhar's mass formula for 
charged stars Q also predicts an infinite mass for 
the extremal case. With a charge to mass ratio a = 
Q/VG^, the Chandrasekhar's mass is [T3 |: 



M, 



ch 



5.83 



(l-a2) 



Ma 



(114) 



This equation reduces to the known Chandrasekhar's 
mass formula when a = (no charge). For the ex- 
tremal case a — > 1, the formula gives a mass tending 
to infinite. 

In the Sub-sec. IVIII Bl we will discuss the tempo- 
ral evolution of neutron stars with extremal electric 
charge. 



2. Binding energy per nucleon 

The gravitational binding energy is given by the 
difference between the total rest mass /i and the total 
gravitational mass m, [s^, p^ : 



B = {fi — m) . 
The binding energy per nucleon is 



B 
A 



[fi — m) (? m 



A 



/ \ 2 



(115) 



(116) 



where m„ is the rest mass of the nucleons, c is the 
speed of light, and A = fi/mn is the total number of 
nucleons on the star. 

The total binding energy of the maximum mass 
model increase with the charge (see Fig. |2K)- 

However, the contrary is true for the binding en- 
ergy per nucleon, i.e. it is lower with higher total 
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charge. For example, the values of B/A for the max- 
imum mass models are: 39.11 Mev per nucleon for a 
neutron star with zero charge; 32.3 Mev per nucleon 
for a star with Q/^/Gfj, — 0.5; 19.63 Mev per nucleon 
for a star with Q/^/Gfi — 0.8; 7.46 Mev per nucleon 
for a star with Qj\fG\x = 0.97, and tending to zero 
as Q approaches the extremal value. For comparison, 
the binding energy per nucleon in the most stable fi- 
nite atomic nucleus is roughly ~ 8 Mev per nucleon 
(see, for example, 0). It could be interesting to 
ask whether a nearly extremal charged star, would 
disintegrate emitting high energy charged nucleus to 
strength its binding energy per nucleon. We will leave 
this question aside by now, since a more realistic EOS 
must be implemented. 

B. Temporal evolution of charged neutron stars 

With the Collapse05v3 it is possible to follow the 
temporal evolution of the stars. We arbitrarily choose 
to study the evolution of three models (8, 13, 21) con- 
structed with the HEOSvl for each of the charge to 
mass values: Q/\/G^ = 0, 0.5, 0.8. Some of these 
models are on the stable branch, while others are on 
the unstable branch (see Fig. [SJ. The region passing 
the first maximum of the mass, where dM/dp < 0, 
corresponds to unstable models at first order in per- 
turbation analysis. We do not study the evolution of 
models passing the second mass maximum (see Fig. 
and tables nnil). 

The results of the simulations are summarized in 
the table Hvl 

It was checked that the stars on the stable branch, 
say for example model 13 with Q/^/Gp. = 0.8, could 
not jump to the unstable branch given a strong initial 
perturbation to the star. The stable stars oscillate 
when perturbed and its (fundamental) frequency of 
oscillation was calculated (see table II V|). 

In the cases where a star oscillates we follow its evo- 
lution over several periods of time to be sure that the 
equilibrium is not a metastable state. For a few mod- 
els, we simulate the oscillating star over roughly one 
minute of physical time (that corresponds to several 
hours of computer time). Over this period of time, 
the amplitude and the speed of the oscillations are re- 
duced due to the numerical viscosity. Of course, the 
simulation can be extended so far as wanted in time. 
Eventually, the velocity will be zero everywhere, and 
the star will rest in perfect hydrostatic equilibrium. 

The models number 8 and 13 are stable and they os- 
cillate when perturbed. We found that the frequency 
of oscillation is higher for lower total charge (see table 

Some of the stable charged stars collapses after a 



sudden discharge of its electric field (see subsection 
below, and table II V|) . 

The charged stars on the unstable branch collapses 
in agreement with the relativistic stellar perturbation 
theory (see for example P^D- 

It is possible to simulate the formation of an ap- 
parent horizon (and a trapped surface) and when this 
happens it is assumed that the star collapsed (see the 
section lVIII B It . In order to keep the accuracy of the 
results the simulations are stopped when an apparent 
horizon forms (see the section IVIip . 

As we will see in Sec. (jVIIIB l|l . for each mass coor- 
dinate there is a value of the coordinate i?, denoted 
as i?+: 

Gm{t, fi) G I ^ Q'^i^i) 
R+{i^N^ -2 + ^Y"^ (^.W (117) 

such that a trapped surface forms at coordinates 
{t,fi,9,(t)) m R{t,fi) < R+{t,fi) (see Sec. IVIIIB1|) . 
In the Eq. (|117ll . m{t, fj,) and Q(^) are the total mass 
and total charge at coordinate /Lt, respectively, G is 
the gravitational constant and c is the speed of light. 

If the perturbed energy is such that the total bind- 
ing energy of an unstable star is greater than its equi- 
librium value, the star will expand first and later col- 
lapse. On the contrary, if the perturbed binding en- 
ergy of the unstable star is lower than its equilibrium 
value, the star will collapse directly to a black hole. 
The perturbed stable stars, will evolve through the 
path that carries it to the corresponding equilibrium 
point on the equilibrium curve (see Fig. |5J. As the 
star contracts or expands its central density changes 
accordingly, and its evolutionary path is an horizontal 
segment in the graph of the total mass versus central 
density, or total binding energy versus central den- 
sity (Figs. EJ). For the collapsing models this segment 
will extend to higher densities until the simulation is 
stopped. 

All the models constructed with the HEOSvl code 
and tested with the Collapse05v3 code have a charge 
Q < VCfis, but for a few models we increase the 
charge to the extremal value Q = VGiIs (conserv- 
ing the total energy). This is easily done on the 
Collapse05v3. The result is that the star exploded, 
resulting in an outward velocity at all the Lagrangian 
points (with the exception at the coordinate origin 
where the boundary condition is maintained: u{fi = 
0) = 0). This is another confirmation of the results ob- 
tained with the HEOSvl code (see subsection lVIII AT]i 
that is not possible to get a finite mass star with ex- 
tremal charge and bounded in a finite spatial region. 

The models number 21 are in the unstable branch 
(see Fig. and we found that all of them collapse 
(see table IIV|l independently of the amount of charge 
(with Q < \/G^s)- However, the time at which the 
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apparent horizon forms is higher for larger total charge 
in the star. In all the cases studied there is no ejected 
matter. 

In the figures it is shown the effect of these pertur- 
bations by means of a series of temporal snapshots. 
The Fig. Q corresponds to the collapsing model 21 
with Q/VG^i = 0.5. In the Fig. ^) it is plotted 
a series of temporal snapshots of the velocity profiles 
for this model. We see that the in-falling matter ac- 
quires a relativistic speed and there are not strong 
shock waves formed. The Fig. |0Jd) shows a series 
of temporal snapshots of the metric coefhcient gu as 
function of the coordinate /x. The metric coefficient 
gtt goes to zero in this case, which is a sufficient con- 
dition for the formation of an apparent horizon (see 
Sec. I VIII B 111 . The Fig. ^ corresponds to the model 
21 with QjyfGpL — 0.5. This is another example of 
an star that collapses without forming strong shock 
waves (see Fig. (SJi)- The results shown in the Fig. 
(|3Jd) also assures that the star collapsed to a charged 
black hole (see Sec. IVIIIBH . 

The Fig. © shows a comparison of the evolution 
of the factor F for two simulations of the model 21 
with charge to mass ratios Q/VG^ — 0.8 (Fig. 
and Q/VGii 0.5 (Fig. Wp)- When the function F 
approachs zero the fluid enters a regime called of "con- 
tinued collapse" , in this case the gradient of pressure 
and of electric charge is no more effective in coun- 
terbalancing the gravitational attraction, the fluid is 
almost in free fall and nothing can stop the collapse. 
From the Eq. (|112a|l we see that the first term inside 
the brackets is nearly zero in this regime. Moreover, 
the figures and ^p) indicate that the factor F 
can acquire negative values after the formation of an 
apparent horizon, in this case the Eq. H112a|l indicate 
that this behaviour reinforce the collapse. 

The Fig. ||7J) corresponds to the oscillating model 
13 with Q/VG/i = 0.5. In this case the amplitude 
of the velocity profiles are bounded, roughly between 
the values ±3000 km sec"^ (Fig. The factor gu 

shows very little variation in this case (Fig. \7}p)- The 
Fig. ((SJ, shows the temporal evolution of five layers 
of the star, which oscillate over several periods. The 
mass enclosed by the layers is indicated in the figure. 



totically flat space-time is the region from where no 
causal or light signals can reach 1+ (the future null 
infinity, see for example [s^ and |2^ for definitions). 
The event horizon is the boundary of the black 
hole region = J~(X+). 

The formation of an apparent horizon is a sufficient, 
although not necessary, condition for the formation 
of a black hole (see Ref. |2^]). In fact, in the case 
an apparent horizon forms the theorems of Hawking 
and Penrose (1970) can be applied to know that the 
outcome will be the formation of a singularity, i.e.: 
the spacetime contains at least one incomplete time- 
like or null geodesies (see the theorems 9.5.3 and 9.5.4 
in ^32], Pags. 239-241, known as "singularity theo- 
rems"). Moreover, the trapped surface (and the ap- 
parent horizon) is contained within a black hole as 
a mathematical proposition asserts (see the proposi- 
tions 12.2.2 and 12.2.3 in Ull, Pags. 309-310, and 
proposition 9.2.1 in P^-S- 311). 

So, we assumed that if an apparent horizon is 
formed in the stellar collapse, the result will be the 
inexorable formation of a black hole. 

In a numeric simulation, in general, is needed an 
algorithm to find out the apparent horizon's loca- 
tion. However, in the present study it is easy to guess 
that the trapped surfaces are spheres. With this in 
mind, we can describe the apparent horizon forma- 
tion and evolution. The vector l'^ = ((ac)~^, 0, 0) 
is tangent to the outgoing null geodesies, and n'^ = 
((ac)^^, — 6^^, 0, 0) is tangent to the ingoing null 
geodesies, in co-moving coordinates. They are also 
orthogonal to the constant {t, fi) surfaces. Following 
Mashhoon and Partovi [23|, we define the quantities 
^I* ~ l^R.i, and $ = n^R^i, which are related by a pos- 
itive multiplicative factor to the expansion factor for 
radially ingoing and outgoing null geodesies, respec- 
tively. Thus, the coordinates of the apparent horizon 
are found solving the equation: 



*(t,Ai) = {ac)-^Rt+b-^R^=Q, 



(118) 



using the Eqs. 146() and (|47|l . and squaring, this is 
equivalent to 



F' = ^ 



(119) 



1. The formation of an apparent horizon 

A trapped surface is a space-like orientable two di- 
mensional compact surface, such that inward and out- 
ward null geodesies normal to it converge. The appar- 
ent horizon is the outer boundary of all the trapped 
surfaces (see [Sll, [2lj, (s^, [H, |2g, and [13, for defi- 
nitions and related theorems) . A black hole in asymp- 



and using Eq. H48|l . this gives: 

2mG GQ^ 
Rc^ d^K^ 

The solution of this equation is: 



(120) 



R±lt,^l) = 



G 



(121) 
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where i?_ is the inner boundary, while is the outer 
boundary of the trapped surfaces. So, the apparent 
horizon is the surface (in Schwarszchild coordinates): 

Sah -.R^ R+, = 6, (j) = (j). 

In order to know the space-time coordinates of the 
apparent horizon in co- moving coordinates, we must 
solve the Eq. H121|l for t and /i. This is done with a 
numeric subroutine that for each time t and for each 
coordinate ^ asks if R{t,fi) = R+{t,fj,) [s^. 

Some of the stars when evolved in time give rise to 
the formation of an apparent horizon, that is the co- 
ordinate R contracted to the value R+ for some t and 
IJ,. The table IIVI summarize the results of the simu- 
lations: it is shown that some of the stars collapsed, 
and it is given the time elapsed from the beginning 
of the simulation and the formation of the apparent 
horizon. 

The Eq. (|118|l is also equivalent to: 



this equation can also be obtained from the equation 
for light rays ds^ = 55] . From this equation we see 
that outgoing light rays have zero expansion in the 
co-moving frame when a = gu ^ |56l |. However, 
from the simulation was obtained that the apparent 
horizon forms before a = 0, this means that radial 
light rays enter the trapped region (at its boundary 
Sah) with d^/dt ^ 0. 

In the Figs. and is shown the evolu- 

tion of the metric coefficient a for model number 21 
with Q/\/Gii = 0.5, and for model number 21 with 
Q/VCfi = 0.8 (this model is one of the entries of the 
table lll|) . As we see, a — > for this models, so they 
collapse. 

The event horizon can not be described in the sim- 
ulations, since it is the result of the complete history 
of the collapsing matter (while the simulation lasts a 
finite amount of time). But we know that the trapped 
region is contained in the black hole region. There ex- 
ists only one situation in which the formation of the 
event horizon can occur in a finite coordinate time: 
when all the matter composing the star collapses in 
a finite time. Only in the special case that the star 
collapse completely, i.e.: if i?(to,/is) < R+{to,^s) for 
some finite time to, the apparent horizon will coin- 
cide with the event horizon of a Reissner-Nordstrom 
space-time: Rbh{M,Q) = R+{to,Hs), where Rbh is 
the Reissner-Nordstrom radius of a black hole of to- 
tal mass-energy M and total charge Q. In all the 
simulations performed the apparent horizon is formed 
for some coordinate /i with < ^ < /is- So, in the 
present set of simulations, the apparent horizon never 
coincides with the event horizon. 



2. Evolution of suddenly discharged stars 

We picked up some stable stars and followed their 
temporal evolution after a sudden discharge of its 
electric field. The scenario is a pure hypothetical 
one, in which the neutron star is leaved in a charged 
metastable state after its formation, or during an ac- 
cretion process onto it. The charged matter could 
suffer a discharge by recombination of charges of op- 
posite sign, or by effect of the high electrical conduc- 
tivity. We emphasize that we are not claiming that 
this is a possible astrophysical scenario, but we want 
to explore the different theoretical possibilities for the 
stellar dynamics. 

The simulation is performed with the worst condi- 
tions for the stellar stability: the electric field is sim- 
ply turned off after the first time-steps of temporal 
evolution of an otherwise stable star. The energy is 
maintained constant in the process, corresponding to 
a conversion of the electromagnetic energy into heat 
or internal energy. This is easily performed with the 
Collapse05v3 code. 

The results of the simulation of the stars suddenly 
discharged are summarized in the table Hvl The mod- 
els discharged are indicated with primed numbers. We 
see that not all the models collapse to form black 
holes: the model number 13 with an initial charge 
to mass ratio Q/VCfi = 0.5 oscillate, while the same 
model with QI\fG[i = 0.8 collapsed after a sudden 
discharge; the model number 8 oscillate if its initial 
charge to mass ratio is Q/\/G^, = 0.5, while the same 
model collapses if QI\fG[i = 0.8. The models num- 
ber 21 always collapse, as they fall on the unstable 
branch, then it is unnecessary to consider them. 



3. The formation of a shell 

In the Reference it is described the formation 
of a shell of higher density formed near the surface of 
the star. Although this effect must happen in Newto- 
nian physics, its evolution in the strong field regime is 
highly non-linear and far from obvious. 

The weight of the star is supported by the gas 
pressure and by the Coulomb repulsion of the mat- 
ter. If the Coulomb repulsion is important respect to 
the gravitational attraction -although not necessarily 
stronger- a shell of matter can form (see for an ex- 
planation). The contrast of density between the shell 
and its interior depends on the energy and charge dis- 
tribution. So, it is important to check if the effect take 
place in the present simulations. 

We found that the shell forms only mildly, or do not 
form, and its late behavior can not be followed clearly. 



17 



From all the numerical experiments performed, it 
is observed that the shell formation effect arise more 
clearly when the initial density profile is flat (like the 
simulations of Ref. consequence, the shell 

formation is not an important physical effect in the 
collapse of charged neutron stars. However, it could 
be important for the collapse of the core of super- 
massive stars 01 . 



IX. FINAL REMARKS 

It is usually assumed in astrophysics that stars 
haven't important internal electric fields. Whether 
an star can have large internal electric fields or a net 
total charge is not yet clear. However, several features 
could be common to the evolution of rotating collaps- 
ing stars, with the angular momentum playing the role 
of the electric charge. So, the present study can shed 
some light on more realistic astrophysical scenarios. 

In the models studied in this paper we assumed, for 
simplicity, that the charge density is proportional to 
the rest mass density. We found that in hydrostatic 
equilibrium the charged stars have a larger mass and 
radius than the uncharged ones. This is, as expected, 
due to the Coulomb repulsion. The mass of the mod- 
els tends to infinity as the charge approaches the ex- 
tremal value Q = VGfis- The hydrostatic equilibrium 
solutions with Q > VCfis gives models with an infi- 
nite mass and radius (see also This means that 
in this particular cases the integrated mass and radius 
diverges within the computer capacity. 

All the models with charge less than the extremal 
{Q < VCfig), have a mass limit and there are unsta- 
ble and stable solutions. We checked the stability of 
the solutions integrating forward in time the models, 
and applying strong perturbations to them. Some of 
the models collapse directly to form black holes, with- 
out ejecting matter. Other models oscillate. For a 
given model, with fixed central density, the frequency 
of oscillation is lower when the charge is higher. The 
frequency of oscillation is weakly dependent on the 
numerical viscosity. The models that collapse are so- 
lutions with dM/dp < 1, while the oscillating models 
have dM/dp > 1. So, the stability of the models agree 
with the predictions of the first order perturbation 
theory [sg. 

It seems that there is a limit for the charge that a 
star can have, which is lower than the extremal case. 
This limit arise because the binding energy per nu- 
cleon of the models with Qj^fGp > 0.97, is lower 
than the binding energy per nucleon on an atomic nu- 
cleus. Then, it is possible that this stars disintegrate 
to reach a more bounded energy state. This point 



deserves further study. 

From a pure theoretical point of view, the issue of 
the collapse of charged fluid spheres, is related to the 
third law of black hole thermodvnamics [l3l , and with 
the cosmic censorship hypothesis jl^j |l9j |. The third 
law of black hole physics states that the temperature 
of a black hole cannot be reduced to zero by a finite 
number of operations. The impossibility of transform- 
ing a black hole into an extremal one, in a finite num- 
ber of steps, is related to the impossibility of getting 
Q = VGhs in some particular experiment 0|, p^ . 
In agreement with this law, we found that it is not 
possible to form extremal black holes from the col- 
lapse of a charged fluid ball. In fact, any charged ball 
with a charge to mass ratio greater or equal than one 
explodes, or its matter spreads. 

The black holes with Q > VGhs represent naked 
singularities, and the impossibility of getting black 
holes with this charge to mass ratio in the present 
simulations are in agreement with the "cosmic cen- 
sorship hypothesis" . 

It must be remarked that once an apparent horizon 
forms the formation of an event horizon is inevitable, 
and as we proved the matching of the interior solu- 
tion with an exterior Reissner- Nordstrom solution, we 
simulated here the formation of a Reissner-Nordstrom 
space-time from a gravitational stellar collapse. 

Although we are not using a realistic equation of 
state we think that the results are of general validity, 
at least from a qualitative point of view. 

Other fields or more exotic physics must be consid- 
ered in order to form extremal black holes. 



X. APPENDIX A 

For completeness we reproduce here the Einstein 
equations in co-movin g co ordinates as derived by Lan- 
dau and Lifshitz (see j23|, Pag. 311). The equations 
of section UTTI were derived from these by making the 
replacements a^c^ = e'^c^, &^ — e^, and — e^. The 
Einstein equations are: 

8^G 1 SttG/ Q2 X 1 /^y2 

-^T^^—[-P+-^]=-e \—+pv 



47ri?4 ) 2 ■ 
1 . . 3 



SttG 



SttG 



6c^ 



-A 



-2/i' — pp' + \p' + v' p 
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TABLE I: Mass, radius, and central density for neutron 
star models with zero charge. 

Model Radius" Mass' Density'^ 
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"Radius at the surface of the star measured in [km] . 
'Total mass-energy (Eq. I45i for each stellar model given in 
solar masses. 

■^Central density of each model measured in [gcm~^]. 

''This is the maximum neutron stars mass, calculated with 
mass steps of 0.0015 Mq. It is also the first maximum indicated 
in the Fig. j^Jj). 

■^This is the second maximum indicated in the Fig. 
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TABLE II: Mass, radius, and central density for neutron 
star models with Q/VG^ = 0.8. 

Model Radius" Mass' Density'^ 
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"Radius at the surface of the star measured in [km] . 
'Total mass-energy (Eq. 1451 for each stellar model given in 
solar masses. 

■^Central density of each model measured in [gcm"^]. 
''This model corresponds to the first mass maximum. 
•^This model corresponds to the second mass maximum. 
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TABLE III: Mass, radius, and central density for neutron 
star models with Q/VG^ = 0.97. 

Model Radius" Mass' Density'^ 
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"Radius at the surface of the star measured in [km] . 
'Total mass-energy (Eq. 1451 for each stellar model given in 
solar masses. 

■^Central density of each models measured in [gcm"^]. 
"^This model corresponds to the first mass maximum. 
•^This model corresponds to the second mass maximum. 
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TABLE IV: Charged neutron star models evolved forward in time. 
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"The primed numbers indicate models that were suddenly dis- 
charged in the simulation. 

''The oscillation frequency of the stable models depends very 
weakly on the numerical viscosity. When the numerical viscos- 
ity coefficient is divided by a factor of two, there is no appre- 
ciable change in the frequency. 

•^The collapse time, is the time elapsed between the formation 
of the apparent horizon and the beginning of the simulation. 

"^Any model with Q/VG/* > 1.0 explodes, i.e., the matter 
scatters to infinity. 
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FIG. 2: Total binding energy in units of Mq(? , for models with different values of the charge to mass ratio QjyfGii 
0, 0.5, 0.8, 0.97 (Fig. a), and mass versus central density for models with different charge to mass ratio Q/VCfx 
0, 0.5, 0.7, 0.8 (Fig. b). 
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FIG. 3: Coefficient grr of the metric for models number 
13 and for three different values of the charge to mass ratio 
Q/VGh = 0, 0.5, 0.8. 
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model 21,Q/(G"» = 0.8 
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FIG. 4: Temporal snapshots of the speed u, for different time-steps (Fig. a), and snapshots of the coefficient gu of the 
metric, for different time-steps (Fig. b), for model 21 and for a charge to mass ratio Q/y/G^ = 0.8. 
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FIG. 5: Temporal snapshots of the speed u, for different time-steps (Fig.a) and snapshots of the coefficient gu of the 
metric, for different time-steps (Fig. b), for model 21 and for a charge to mass ratio Q/VOfi = 0.5. 
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FIG. 6: Temporal snapshots of the factor F, for different time-steps, for model 21 and for a charge to mass ratio 
Q/\/G^, = 0.8 (Fig. a), and temporal snapshots of the factor F, for different time-steps, for model 21 and for a charge 
to mass ratio Q/VG^l = 0.5 (Fig. b). 




FIG. 7: Snapshots of the speed u (Fig. a) and snapshots of the coefficient gtt of the metric (Fig. b), for different 
time-steps over half period of oscillation, for model 13 and for a charge to mass ratio Q/^/G^ = 0.5. 
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model 8, Q / (G°= n) = 0.8 
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FIG. 8: Temporal evolution of several layers of an oscillating star, for model 8 and with charge to mass ratio Q / VOfx = 0.8, 
over several oscillation periods. Each curve is labeled by the quantity of rest mass enclosed, or equivalently by its 
Lagrangian coordinate n, given in units of the total rest mass /Xs 



